The fast sampling algorithm for Lie- Trotter products 
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A fast algorithm for path sampling in path integral Monte Carlo simulations is proposed. The 
algorithm utilizes the Levy-Ciesielski implementation of Lie- Trotter products to achieve a math- 
ematically proven computational cost of nlog 2 (n) with the number of time slices n, despite the 
fact that each path variable is updated separately, for reasons of optimality. In this respect, we 
demonstrate that updating a group of random variables simultaneously results in loss of efficiency. 
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Path sampling in path integral Monte Carlo simula- 
tions becomes more difficult in the limit of a large num- 
ber of path variables or time slices, not only because 
of the build-up of correlation among the variables that 
are sampled, but also because of the shear increase in 
the number of random variables. Over the time, various 
approaches have been attempted in order to overcome 
the slowing down of the simulation. For direct sampling 
of Lie- Trotter products, among the most successful are 
the staging methodQ, the threading algorithm 0, the 
bisection method U and the multigrid technique 0. 
The so-called "normal mode" and Fourier formulations of 
ath integrals have also been shown to improve sampling 
"l 0] • The last techniques are part of a larger class 
of path integral methods, class that is called the random 
series implementation |8(]. To give a few examples, the 
computational effort in most random series approaches 
scales as n 2 with the number of path variables, whereas 
the bisection method may reduce the effort down to n 1A 



The random series approach, at least in the primi- 
tive form, is not particularly efficient However, it 
reveals the incredibly large variety of possible different 
path integral formulations, while also being suggestive of 
more optimal approaches 9] . At the same time, it shows 
the strong connection that exists between these differ- 
ent forms, connection that is realized by means of cer- 
tain orthogonal transformations. For Lie- Trotter prod- 
ucts, Predescu and Doll have shown that there is an 
infinity of possible normal mode representations, which 
can be obtained one from the other by certain orthogo- 
nal transformations. One such transformation leads to 
the Levy-Ciesielski representation, which has very spe- 
cial properties when it comes to numerical implemen- 
tation. For instance, it allows for fast computation of 
paths, with a scaling of nlog 2 (n) for an entire path, as 
opposed to n , for other representations. The only re- 
striction is that the number of time slices must be of the 
form n — 2 k , a condition reminiscent of the fast Fourier 
transform. The algorithm we propose for the sampling 
of Lie- Trotter products utilizes the Levy-Ciesielski repre- 



sentation to achieve a similar goal: sampling in n\og 2 (n) 
operations of entire paths, while updating each path vari- 
able individually. Since any algorithm that computes en- 
tire paths every time has a computational effort propor- 
tional to n, we see that the proposed algorithm is highly 
efficient. Since it is not necessary to update all path vari- 
ables all the time, the proposed algorithm can be made 
even more efficient by using it in conjunction with the 
multilevel Monte Carlo sampling technique jj]. 

We commence by casting an arbitrary Lie- Trotter 
product in the Levy-Ciesielski form. In order to do so, we 
first enumerate the basic properties of the Levy-Ciesielski 
series representation of the Brownian bridge. For more 
details, the reader is advised to consult Refs. 0, fill Il2t 
For k = 1,2,... and j = 1,2,..., 2 fc ~ 1 , the Schauder 
functions Fkj(u) are generated by translations and di- 
latations of the function 



u, ue (0,1/2], 
F 1A {u) = { 1-u, we (1/2,1), 
0, elsewhere. 



More precisely, we have 



F fcJ (u) = 2-( fc - 1 )/ 2 F 1;1 (2 



k-X„ 



■J + l), 



(1) 



(2) 



for k > 1 and 1 <j < 2 fc ~ 1 . 

If we multiply them by 2 _ ( fe -i)/2^ the Schauder func- 
tions make up a pyramidal structure organized in layers 
indexed by k, as shown in Fig.^ The supports (the sets 
on which the functions do not vanish) of the Schauder 
functions are the open intervals of the form (ukj-i, Uk.j), 
for 1 < j < 2 fc_1 , where Uk,j = j2~( k ~ 1 \ The supports 
are disjoint for functions corresponding to the same layer 
k. Because of this property, we have the equality 

2fc -i 

a KjFk,j(u) = a fc [ 2 fc-i tl ] +1 i ;i fc [2 fc-i u ] + i(u), (3) 

for any sequence of numbers ak,i, dfc,2, ■ ■ ■ , Ofc^- 1 ■ Here, 
[x] denotes the largest integer smaller or equal to 
x, whereas for u = 1, the quantities 2 h - t +i an d 
Fk 2 fc - 1 +i(l) are defined to be equal to 0. 
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FIG. 1: A plot of the renormalized Schauder functions for 
the layers k — 1,2, and 3, showing the pyramidal structure. 



Let {dk,j',k — 1,2, ... ;j — 1,2,..., 2 fe ~ 1 } be an infinite 
sequence of independent identically distributed standard 
normal variables. From Eq. © and the Leyy-Ciesielski 
construction of the Brownian bridge |10L llll Il2| , we have 
that 



B< 



z , afc,[2 fc 
fc=i 



1 u] + l- t k,[2 k - 1 u]+l 



(4) 



In words, the right-hand side random series is equal in 
distribution to a standard Brownian bridge. 

Let n — 2 k — 1 be a fixed number and consider the 
equidistant points Uj = j2~ k , with 1 < j < 2 k — 1 = n 
(these points were denoted before by Uk+i,j, but we shall 
drop the index fc + 1 to avoid cluttering the formulas). 
Because the Schauder functions Fi : i(u) vanish at these 
points for all levels I > k + 1, it follows that the random 
sums 



i=i 



a l,[2'- 1 u j ]+lPl,[2'- 1 u j ]+l( u j) 



for j = 1,2, ... ,n have joint distribution equal to the 
joint distribution of the random variables B® L . . By the 



definition of the Brownian bridge, the joint distribution 
of the latter variables is given by the formula 



1 



(0,xi)p U2—U1 {%1 1 ^2) ' ' ' Pl — u n (*n,0), (5) 



Pi(0,0) J 

with p u (x, x') defined by 

p u (x, x') = (2ttu)- 1/2 exp [-(x' - x) 2 /{2u)] . (6) 

Using this observation, the notation x r (u) = x+(x' — x)u, 
and the definition a 2 — h 2 /3/m , one easily proves that 
the joint distribution of the random variables 



.(Uj) +0"^ a/,[2 i - 1 « J ] + l^/,[2 i - 1 M J ] + l( l 



1=1 



for j = 1, 2, . . . , n is given be the formula 

(X1,X2) 

■ ■ ■ Pa 2 (l~u n ){Xn,x') /p a 2 (x, x') 



(7) 



(8) 



Because the density matrix of a free particle is strictly 
positive, any short-time approximation po(x,x';f3) can 
be put in the product form 



p Q (x, x';(3)= Pf p (x, x'\ (3)r Q (x, x'; (3). 



(9) 



Letting xq = x, x n +i — x', uq — 0, and u n +i — 1, the 
n-th order Lie- Trotter product obtained from the short- 
time approximation considered above takes the form 

~ n 

p n (x,x';(3) = / Y\_p a 2 {u ._ u . +l) (xi,x i+1 ) 

n 

x Y[r (x :j ,x : j +1 ;f3/2 k )dx 1 ■■■dx n . (10) 

From the last equation and the fact that the distribu- 
tion given by Eq. (JSJ) is the distribution of the random 
variables appearing in Eq. Q, we readily obtain the fol- 
lowing Levy-Ciesielski form of Lie- Trotter products for 
n + 1 = 2 k time slices 



p n (x,x';(3) = p fp (x,x';P) / da 



[ f da k , 2k - 1 (2TT)- n / 2 l[l[exp(-a 2 l /2 



X il r » 

3=0 



1=1 i=l 



X r (Uj) + a y^a;J 2 i-l Mi | + 1 F;j2'- 1 M,1 + l("j)7 



1=1 



c r (u J+ i) +a ^a; ! [2 ! - 1 « J+ i]+i^,[2 ! - 1 « J+ i]+i(' i j+i);/3/2 fc 



1=1 



(11) 



3 



This formula has been first introduced in Ref. albeit 
for some specialized short-time approximations. 

We have already mentioned that one of the advantages 
of the Levy-Ciesielski form for Lie- Trotter products is 
that it enables fast computation of paths |lfj|. while main- 
taining the random series appearance of the final expres- 
sion. This is so because for each Uj, one needs to perform 
a number of fc = log 2 (n + 1) operations in order to com- 
pute the coordinate Uj of the path. This translates into a 
scaling of (n+ 1) log 2 (n + 1) operations for a whole path. 

As announced in the beginning of the letter, a second 



No two factors defined by the curly brackets contain a 
same variable aij. Indeed, given i G {1, 2, . . . 2 i_1 }, only 
the function Fi^(u) may be non-zero for the values Uj 
with (i - l)2 k -' l+1 <j< i2 k - l+1 . Thus, each factor de- 
fined by the curly brackets appearing in Eq. (|12|l contains 
one and only one variable ai t %. Therefore, the variables 
ai t i are and should be treated as independent during the 
Monte Carlo simulation. Each proposal aij — * a[ i must 
be tested separately using the weight given by the appro- 
priate factor and accepted or rejected according to the 
Metropolis-Hastings rule. 

Let us analyze the efficiency of the algorithm. First, 
there are fc = log 2 (n + 1) layers. For each layer, one 
evaluates the function tq(x, x'\ (3) exactly n+1 times, in 
order to update all variables from the layer, individually. 
Thus, the computational effort to update all variables in- 
dividually is proportional to (n+1) log 2 (n+l). Since any 
algorithm must have a scaling of at least n + 1 [this is 
the computational effort necessary to evaluate the distri- 
bution given by Eq. (|llf) for any update attempt] , we see 
that the technique is extremely efficient, especially given 
that all variables are moved separately. Of the variables 
ai t i, the most difficult to sample is ai t i, because its distri- 
bution stretches over a larger region [the distance covered 
by a variable aij is of the order o-2~( l ~ 1 ^ 2 , for low tem- 
peratures]. Therefore, the overall computational effort to 
perform an accurate sampling of all variables also scales 
as (n + 1) log 2 (n + 1). 

The reader may ask why we have insisted on up- 
dating each path variable individually. The answer is 
that there is a loss of efficiency if we try to update 
more than one path variable at a time. We shall prove 



advantage of the Levy-Ciesielski form is that it allows 
for fast Monte Carlo sampling of the distribution given 
by the right-hand side of Eq. (f 1 1|) . The algorithm is as 
follows. A trial move is proposed for all 2 l ~ x variables 
ai t i that correspond to a single level I. However, the ac- 
ceptance/rejection decision is taken individually for each 
path variable, because these variables are statistically in- 
dependent. To see this, first observe that the part of the 
product distribution in Eq. l]l l|l that only involves the 
path variables for the layer / factorizes as 



(12) 



this assertion in the remainder of the letter. Let's as- 
sume we are given a finite collection X\ , A 2 , . . . , X n 
of independent identically distributed random vectors 
(i.i.d.r.v's), taking values in some space K d . Let p(x), 
with x £ M. d , be the normalized distribution of any of 
the random vectors Xi. By independence, the overall 
distribution is given by the product p(xi)p(x 2 ) . . . p(x„), 
which is a distribution on the space M. dn . Assume we 
attempt to update all variables at once, using a trial dis- 
tribution T(y 1 |x 1 )T(y 2 |x 2 ) . ..T(y n |x„). The move to 
(yi,y 2 , . . . ,y„) is then accepted with probability 

and rejected with the remaining probability. The average 
acceptance probability is given by the formula 

Ac(n) = / dxidyi--- / dx„dy„p(xi)T(yi|xi) 

jR 2d JWL 2d 

■ ■ ■ p(x„)T(y„|x„) min { 1, f[ ^f^] ) .(14) 

In these conditions, we have the following theorem that 
guarantees that simultaneous sampling of i.i.d.r.v's is in- 
efficient. 

Theorem 1 (Bad sampling of i.i.d.r.v's) Except for 
the ideal case T(y|x) = p{y) whenever T(y|x) ^ 0, there 
is a strictly positive constant H such that 

Ac(n) - e- Hn . (15) 



2& — i i i2 k ~ [Jrl — 1 

n \ ^r i/2 e-< j2 n ^ 

j"=(j-l)2 fc - ! + 1 



i=l 



.(uj) + a y^a;j2'-i^i+i-Pi,r2 i - 1 « j i+i(' 1 



i=i 



x r (u j+1 ) + a^a lA2 i-i Uj+l]+1 F tA2 i-i Uj+l]+1 (u :j+1 ); (3/2 k 
I 



i=i 



4 



This constant is given by the relative Shannon entropy 

>(y)T(x|y) 



H = - 



p(x)T(y|x)log 



p(x)T(y|x) 



dxdy. (16) 



Proof of the theorem. For convenience, we let E denote 
the expected value against the distribution p(x)T(y|x), 
i.e., for some arbitrary function /(x, y), 



E[f(X,Y)]= / p(x)T(y|x)/(x,y)dxdy. 
Then Eq. (|14f> can be written as 

p<Xi)T(Xi\Yi) \ 



Ac(ra) = Ei 



(17) 



Now, consider the identity 

\ p(Xi)T(Yi\Xi) 
expf-J-i^Sog 



n 



pjY^TjXtlY) 

P (Xi)T(y i |Xi) 



(18) 



The expression inside the curly brackets is a "time" aver- 
age of independent identically distributed random vari- 
ables. By the law of large numbers, this time average 
will converge to a constant function, the value of which 
is the "space" average 



H 



lim 

n — >oo 77, 



1 " 

- Y] log 



E<Mog 



pjY^TjXtlY) 
p{Xi)T(Yi\Xi] 

p(Y)T{X\Y) 



P {X)T(Y\X) 



(19) 



Remembering the definition of E, we see that the right- 
hand side of the previous equation is the Shannon en- 
tropy of the probability measure p(x)T(y|x) relative to 
the measure p(y)T(x|y). It is always non-negative and, 
in fact, it is (strictly) positive except for the case 



p(y)T(x| y ) = p(x)T(y|x). 



(20) 



To prove the last assertions, notice that — log(x) is a 
strictly convex function and remember Jensen's inequal- 
ity, which in this case says 

E {- log[/(X, Y))}>- log {E [f(X, Y)}} , 

with equality if and only if f{X, Y) is constant. Then, 

>(y)T(x|y)" 



H = - 
log 



> 



p(x)T(y|x)lo; 
/ p(x)T(y|x)^ 

JR™ P(X 



p(x)T(y|x) 
Ky)T(x|y) 



dxdy 



P(x)T(y|x) 



dxdy 



= 0. 



Eq. (0 follows from Eqs. (EJ, 0, and JEjJ, to- 
gether with the observation that min{l, exp(— Hn)} = 



exp(— Hn), since H > 0. Eq. i|20|) follows from the sec- 
ond part of Jensen's inequality. It is readily seen to be 
equivalent to the statement p(y) — T(y|x), whenever 
T(y|x) 7^ 0. The proof of the theorem is concluded. 

We have therefore demonstrated that there is a loss 
of efficiency if simultaneous updating of path variables is 
attempted, even for uncorrelated variables. If simultane- 
ous sampling is performed, one must modify the proposal 
T(y|x), so that the Shannon entropy decreases with the 
dimensionality at a rate faster than 0(l/n). Straight- 
forward calculations show that, if the proposal T(y|x) 
is uniform in a d-dimensional hypercube centered about 
the previous position, the maximal displacements must 
be decreased at a rate faster than or equal to 0(n -1 / 2 ), 
in order to prevent the severe degradation of the qual- 
ity of the simulation predicted by Th. ^ For arbitrary 
random series, there is little choice: either we decrease 
the maximal displacements or update each path variable 
individually, at a total cost proportional to n 2 . 

However, for the Levy-Ciesielski representation, the 
random variables can be updated individually in 
n log 2 (n) operations, where n is the number of time slices. 
For this reason, as well as for the property of fast compu- 
tation of paths, we believe that the sampling algorithm 
we have presented will prove to be a valuable tool for 
all path integral simulations that are implemented via 
Lie- Trotter products. Again, the main property of the 
Levy-Ciesielski representation that has enabled the de- 
velopment of this algorithm is the fact that the path 
variables corresponding to a same layer are statistically 
independent. 
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